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Abstract 

We show analytically that the inverse kinetic inductance L^ 1 of an over- 
damped junction array at low frequencies is proportional to the admittance 
of an inhomogeneous equivalent impedance network. The ij bond in this 
equivalent network has an inverse inductance Jij cos(9® — 6® — Aij), where 
Jij is the Josephson coupling energy of the ij th bond, 9® is the ground-state 
phase of the grain i, and Aij is the usual magnetic phase factor. We use 
this theorem to calculate L _1 for square arrays as large as 180 x 180. The 
calculated L^ 1 is in very good agreement with the low-temperature limit of 
the helicity modulus 7 calculated by conventional equilibrium Monte Carlo 
techniques. However, the finite temperature structure of 7, as a function of 
magnetic field, is sharper than the zero-temperature -L -1 , which shows sur- 
prisingly weak structure. In triangular arrays, the equilibrium calculation of 
7 yields a series of peaks at frustrations / = ^(1 — 1 /N), where N is an integer 
> 2, consistent with experiment. 

PACS numbers: 74.50.+r, 74.60.Ge, 74.60.Jg, 74.70.Mq 
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I. INTRODUCTION. 



Superconducting arrays have been a subject of continuing interest for a number of years 
|Ij]. Such arrays typically consist of a large collection of superconducting grains, of micron 
or submicron dimensions, arranged in a periodic structure. They are embedded in a non- 
superconducting host which mediates a Josephson or proximity effect coupling between the 
grains. The coupling generally produces a transition to a coherent state at a temperature 
T C (B) which is well below the superconducting transition temperature T c0 of the individual 
grains, and which strongly depends on the applied magnetic field B. There has been exten- 
sive research on the thermodynamic properties of these arrays, especially as a function of 
B . More recently, a number of workers have extended these studies to the dynamical 
response of such arrays, especially their a. c. IV characteristics [H and vortex motions in 



the array [10 



In theoretical studies, the thermodynamic response of superconducting arrays is com- 
monly characterized by the magnetic field and temperature dependent helicity modulus, 
"f(B, T). 7 is defined as the stiffness of the array to a twist in the phase of the supercon- 
ducting order parameter It is thus closely analogous to the spin wave stiffness in a 
magnetic material. 7 can readily be expressed as an equilibrium quantity and evaluated 
as a thermodynamic average in the canonical ensemble, e. g., by Monte Carlo simulation. 



On the other hand, the measurement of 7 is typically a dynamical problem [11,12] . Such 
a measurement is based on the identification of 7 with the inverse kinetic inductance L^ 1 , 
where L is the effective inductance of the array. L can be measured at temperatures below 
T C (B) by applying an appropriate a. c. current and measuring the out-of-phase part of the a. 
c. voltage response. Thus, in order to compare theory and experiment, one usually equates 
a thermodynamic quantity to a dynamical one. 

In this paper, we describe a straightforward dynamical procedure for calculating 7 in a 
superconducting array at temperature T = 0. The method starts from the coupled Josephson 
equations for the array in the limit of low frequency and small amplitude. By carrying out 
an appropriate Taylor expansion in this limit, we show that 7 is equivalent to the inverse 
inductance of a certain inhomogeneous impedance network. The impedances in this network 
are determined by the ground-state phases arrangement in the array. Once the ground state 
is known, the effective impedance of the network can be calculated by standard numerical, 
or even analytical, techniques. 

To show that this method is practical, we apply it to a square Josephson-coupled network 
subjected to an applied transverse magnetic field. The resulting j(T = 0) has surprisingly 
weak structure as a function of magnetic field. Remarkably, the structure appears most 
conspicuously at finite temperature, because T C (B) is an extremely sensitive function of 
B. To verify the consistency of this picture, we calculate 7 by the standard equilibrium 
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approach. We show that, for all values of B considered, it approaches the T = value 
obtained by dynamical techniques. 

We have also used the equilibrium approach to calculate ^y(B,T) for a triangular lattice 
as a function of magnetic field. The dynamical approach is still applicable, in principle, 
to triangular lattices. However, we have not used it in this case, because our algorithm 
for computing impedances does not always converge for triangular lattices. Our calculated 
equilibrium j(B,T) has a very rich structure, as a function of B, which agrees in some detail 
with that found in recent dynamical measurements by Theron et al This agreement 
confirms, once again, that the static and dynamical quantities are basically equivalent, and 
also shows that the array is well described by the standard frustrated XY model which 
underlies both the static and dynamic theories. 

We turn now to the body of the paper. Section II describes the formalism underlying our 
calculations. In Section IIA, we prove the equivalence of j(T = 0) to the effective impedance 
of an inhomogeneous impedance network. Section IIB reviews the standard equilibrium 
method for calculating 7. Section III describes our numerical method for calculating both 
the effective impedance and the quantities needed for an equilibrium evaluation of j(T). Our 
results for both square and triangular lattices are presented in Section IV. A brief discussion 
follows in Section V. 



II. FORMALISM. 



Our basic geometry is shown in Fig. 1. It consists of a periodic network of superconduct- 
ing grains, each of which is Josephson coupled to its nearest neighbors. We can calculate 
the inverse kinetic inductance of this network by two different methods. In the first 
method, we use an identity, proven below, which shows that L^ 1 is equivalent, within a 
proportionality constant, to the admittance of an inhomogeneous equivalent network. This 
identity permits calculation of L~ l , or equivalently the helicity modulus 7, at temperature 
T = 0. The second method is particularly appropriate at finite temperatures. It uses the 
fact that 7 is an equilibrium thermodynamic quantity which can be calculated by standard 
Monte Carlo techniques. We review this type of calculation in Section IIB below. 



A. Helicity Modulus at T = 0: Dynamical Evaluation. 

We consider an overdamped resistively-shunted Josephson junction (RSJ) array in the 
presence of an external magnetic field B applied perpendicular to the array. The dynamical 
equations for the array at zero temperature take the form 

Iij = + h-ij sin(6»i - 9j - A i:j ), (1) 
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^ hj = h;ext- (3) 
i 

Here iy is the current from grain % to grain j, Vij = Vi — Vj is the voltage difference between 
grains i and j, J c; jj and are the critical current and shunt resistance of the ij th junction, 
9i is the phase of the order parameter on the i th grain, and Ii- ex t is the external current fed 
into the i th grain. is a magnetic phase factor defined by the relation 



-/'A-dl, (4) 

Ji 



where A is the vector potential for the external magnetic field, and <E> — hc/2e is the flux 
quantum. Eq. (1) simply expresses the total current of each junction as the sum of a normal 
current through a shunt resistance and a Josephson supercurrent. Eq. (2) is the Josephson 
relation connecting the voltage drop across a junction to the phase difference. Eq. (3) is just 
the Kirchhoff 's Law, expressing the current conservation at the i th grain. In the calculations 
below, we will assume that the vector potential is that of the externally applied magnetic 
field. This is equivalent to assuming that the intergranular Josephson coupling is so weak 
that the screening magnetic fields can be neglected. 

Let us assume that in the absence of an external current, the i th phase has a value 
of Q®. We assume that the collection of phases {9®} represents a local energy minimum. 
One possible state of this kind is the true ground state of the array, but there may be other 
metastable states which also satisfy this condition. We now consider an a. c. external current 
h-,ext =Re [ijo ex p( — ioot)], where Iio is an appropriate amplitude. If the amplitude is not too 
large, the induced voltages and phase changes will also be small, and we may assume that 
they also have frequency uo. Hence, we write the phase on the i th grain as 

We can now use these assumptions to linearize the equations of motion about the state {9®}. 
Combining Eqs. (1) - (3), we can write 



E (h-, i3 sin(0° - 9° - A tJ + 59, - 59,) + K \ R ^ , - n^t- 



Expanding the sine function in a Taylor series to first order in 59{ — 59 j, and using the fact 
that in the absence of external current, the state {9®} satisfies 



E/^sin(^-^-^,) = 0, (5) 

i 

we obtain the simple relation 
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E 



(60i - sej)!^ cos(e9 - 0j - Ay) + 



h(56i - 59 ■ 



2ei? 



\ext- 



(6) 



Finally, we use the assumption that 50, is varying sinusoidally in time with frequency uj to 
rewrite this equation as 



E 



2ze 

-J c; ,cos( 



- 0? - A«) + 



1 



R 



>j J 



Vi,- 



Ii;exti 



(7) 



2e v 



(8) 



where we have used the relation 

We can interpret equations (7) in a very simple way. Namely, they are identical to those 
for an inhomogeneous impedance network responding to the driving a. c. current. The ij th 
admittance is the sum of two terms in parallel: a conductance g^- = 1/Rij, and a purely 
imaginary inductive element [2ie/ (hu)]I c; ij cos(0° — 6® — Aij). The elements of this network 
are all known if the unperturbed metastable state {0°} is known at the fields of interest. 
Given the elements of the network, the effective network admittance can be obtained by 
standard numerical techniques which have been developed for finding the admittances of 
inhomogeneous impedance networks fll3|Jbl . 



The results are particularly simple in the limit of low frequency. In this case, the shunt 
resistances are irrelevant, and the response of the network is purely inductive. The effective 
inductance of the network is then often known as the kinetic inductance . The present work 
permits this kinetic inductance to be computed dynamically, at least at T = 0, by solving 
for the impedance of an effective impedance network. A numerical technique for carrying 
out this calculation is summarized below. 



B. Helicity Modulus at Finite Temperature: Equilibrium Evaluation. 

The kinetic inductance of an overdamped Josephson network can also be calculated as an 
equilibrium quantity, by connecting it to the so-called helicity modulus of the network. This 
connection has been known for more than ten years |§, but for completeness we summarize 
this relation in the following paragraphs. 

The equilibrium Hamiltonian of the Josephson network is given by 

H = - J2 Eij cos(0, - 9j - A,^) (9) 

ij 

where = hl c ,ij/ (2e) is the coupling energy between grains i and j, and is given by Eq. 
(4). In our calculations, we assume that is the same for each junction. The equilibrium 
thermodynamic properties are given by the free energy 
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T= -k B T In Q N (T) 



(10) 



where the TV-particle partition function Qn is given by the canonical integral 

Qn(T)= •••/ IUdOiexp(-H/k B T). 
Jo Jo 

Similarly, the canonical average of any desired function 0(9±, • • ■ ,&n) is given by 



< O >-- 



Q 



N 



• J UidOiO^, ■■■,9 N ) exp(-H/k B T). 



(11) 



(12) 



The helicity modulus 7 Q/ g (or equivalently, the superfluid density) is defined as the free 
energy cost of imposing a twist in the phase at the boundaries of the sample. The principal 
elements of 7^,3 can be thought of as spin-wave stiffness constants appropriate to "phase 
waves" in this weakly coupled system. Rather than imposing a twisted boundary condition 
and calculating the resulting increase in free energy, however, it is more convenient to use 
periodic boundary conditions and calculate 7 Q/ 3 as 



7a/3 = 



' d 2 T \ 
dA' a dA'J Al=Q 



(13) 



Here A' represents an added uniform vector potential (in addition to that which produces 
the applied magnetic field) which is included in the Hamiltonian in order to produce a 
twist. The various second derivatives in Eq. (13) are readily computed for an ordered or a 
disordered sample, with the result for, e. g., j xx : 



EyXy sin(0j - 6j - A^) 




<ij> 



EijXy sin(0j - 9j - A^) 

<ij> 1 



(14) 



where x^ = Xj — Xi and Xi is the x coordinate of the center of grain i. Similar expressions 
hold for the other components of 7. The above expressions are valid for both square and 
triangular lattices. 

We now connect 7 a/ 3 to the inductive response of the array. To make this connection, 
note that the first derivative of the free energy is the local current density J, i. e., 



dA~ 



= -J a /c. 



a / A'=Q 



Hence, 
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1 (dJ a \ 

""--cm)** 

If A' is produced by an a. c. electric field E of frequency u, then A' = —icE/uo. If we intro- 
duce the conductivity a a /3 by a a p = (dJ a /dEp) g_ , then it follows that 7 a /3 = —iua a i3/c 2 . At 
very low frequencies, a a p is of the purely inductive form a a p = iL~^/co, where = c 2r y a p. 
Thus 7 aj g is indeed proportional to low-frequency inductive response L~l of the array 

In general, 7 Q/ g is a tensor, not a scalar. In square or triangular arrays, however, the 
fourfold or sixfold rotational symmetry implies that 7^ is, in fact, just a multiple of the 
unit tensor, i. e. j a p(T) = j(T)8 a p. 

III. METHOD OF CALCULATION. 

In all the calculations in this paper, we assume that the array is homogeneous and 
ordered. That is, we assume that the coefficients I c -ij and Rij are independent of position 
and equal to I c and R respectively. It follows that the coupling energies are also position- 
independent and equal to a constant which we denote simply as Ej = hl c /(2e). 



A. Effective Impedance. 

For the impedance calculation, we inject a uniform a. c. current I ext into each grain at 
one boundary and extract the same current from each grain in the opposite boundary. To 
obtain the kinetic inductance by the method of Section IIA, we need to compute the effective 
impedance of a certain equivalent network. The inverse inductance of the ij th bond in this 
network is l c cos(0° — 8® — Aij), where 9® denotes the ground-state phase of the i th grain. 
Thus, the problem divides into two parts: (i) finding the ground-state configuration {8®} in 
the absence of current; and (ii), given that configuration, computing the effective impedance 
of the network. In practice, we cannot be assured that we have obtained the true ground 
state; in general, we will obtain a metastable configuration which may be close to the ground 
state. 

To obtain a metastable state near the ground state phase configuration, we carry out 
a standard Monte Carlo annealing [15| at fixed magnetic field, starting from a high tem- 



perature and cooling gradually. The purpose of this annealing is to avoid, or at least to 
minimize, the trapping into metastable states which is a notorious consequence of deal- 
ing with frustrated systems, such as these arrays. In this part of the calculation, periodic 
boundary conditions are used in all directions. To check that we are indeed dealing with 
a near-ground-state configuration, we generally carry out more than 50 anneals at a given 
temperature. The effective impedance is computed for the state with the lowest energy. 
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To calculate the effective impedance of the network, we need to solve the set of sparse 
linear equations (7). In practice, at the low-frequency limit, the shunt resistance contribution 
can be discarded, and the trivial common factor of i/uj divided out. For both square and 
triangular arrays, we then use the Gauss-Seidel iterative method to solve for the impedance 
of networks as large as 180 x 180. This iterative procedure always converges in our square 
arrays calculations. In our triangular array calculations, the iterative approach sometimes 
fails to converge, possibly because the ground state factors cos(6>° — #j — A^) can be of both 
signs for triangular arrays. Thus, we will present results for triangular arrays using only the 
equilibrium method described above. 



B. Helicity Modulus at Finite Temperature. 

To calculate the helicity modulus 7 a/ 3 at finite temperature, we use the equilibrium 



Monte Carlo method within the standard Metropolis algorithm |L6[ for both square and 
triangular arrays with periodic boundary conditions. In order to minimize the effects of 
metastable states, we cooled our system from high temperatures at each magnetic field, 
typically starting from LhEj/ks and cooling to 0.5Ej/kB in steps of O.lEj/kB, followed by 
an additional cooling down to 0.2Ej/kB in steps of 0.05Ej/kB- Following this cooling, the 
required averages were obtained by making 60000 passes through the entire lattices, with 
the first 10000 passes discarded. We carry out the thermal average by including only every 
tenth pass through the entire array, in order to minimize the strong statistical correlation 
between successive configurations. The reported results are averages of 7 = {pf xx + j yy )/2 
and are generally produced by averaging over four to eight independent runs. 



IV. RESULTS. 
A. Square Arrays. 

Fig. 2 shows the helicity modulus 7 plotted as a function of frustration values / at four 
different temperatures in square arrays. The array sizes are either 48 x 48 or 50 x 50 in 
the finite-temperature calculations, which are carried out by the equilibrium technique of 
Section IIB. In the T = calculations, we use arrays from 120 x 120 to 180 x 180 and 
the analog technique of Section IIA. The frustration / = $/$0) where <3? is the flux per 
plaquette. Note that as temperature decreases, these Monte Carlo calculations converge 
remarkably well to the T = impedance results. This supports the validity of our method 
for calculating the kinetic inductance at temperature T = 0. 

For / = and / = 1/2, the melting temperature is known to equal approximately 
0.95Ej/k B and 0A5Ej/k B respectively MM- For / — p/q, with p and q mutually prime 
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integers, as q increases, the melting temperature, estimated as the point where 7 goes to 
zero, decreases very quickly, consistent with previous work ||. At / = 5/12, for example, 
the melting temperature is probably less than O.lEj/kg (the vanishing of 7 at this value of 
/ is considerably broadened by finite size effects). 

A remarkable feature of the results shown in Fig. 2 is that the plot of 7(/) becomes 
sharper as the temperature is increased. This sharpening with increasing temperature is 
seemingly contrary to intuition, but in fact is easily explained. As noted above, the melting 
temperature of the vortex lattice at a frustration / = p/q decreases with increasing q 
(assuming p and q mutually prime integers). Thus, at extremely low temperatures, the 
vortex lattice is stable for a wider range of / values, but at slightly higher temperature, it 
has melted at all except a few values of / corresponding to relatively small values of q. The 
corresponding j(f) shows sharp peaks only at a few such values of /. (As seen in Fig. 2, the 
sharpest features at a temperature of O.lEj/kB occur at / = 1/5, 1/3,2/5, and 1/2.) This 
explains why the zero-temperature 7(/) is surprisingly featureless: for all values of /, the 
vortex lattice is apparently frozen at T = 0. 

It should be emphasized that the T = calculations are carried out using "ground 
states" obtained by Monte Carlo annealing from a finite temperature. As pointed out in 
the previous section, such annealing is not guaranteed to produce the true ground state, 
but only a metastable state close to the ground state. We have checked for metastability by 
calculating 7 at T = using several such annealed states. The resulting 7's are generally 
quite close to that obtained in the state of lowest energy. The 7's plotted in Fig. 2 always 
correspond to this lowest energy state. 

B. Triangular Arrays. 

Fig. 3 shows the calculated helicity modulus 7, plotted as a function of the frustration 
value / in a triangular array at four different temperatures. As in a square array, when 
f = p/q, where p and q are mutually prime integers, the melting temperature tends to 
decrease strikingly with increasing q. Indeed, when q is sufficiently large (typically, when 
q > 8), the melting temperature of the array is less than O.lEj /ks)- As in the square array, 
there are a series of sharp peaks in ^(f, T) plotted as a function of / at fixed T. Some of this 
structure has been previously studied [|J, but the present work represents a considerably 
more detailed investigation. We have found clearly distinguishable peaks in 7 at frustration 
values / = 2'3'4'l'f'l'""' m deed, we have obtained peaks at every value of / between 
1/3 and 1/2 satisfying the relation / = |(1 — 1/N), with N as large as 8. All such peaks 
have been observed in recent experiments by Theron et al [Q. Since we are able to obtain 
such peaks here, we conclude that the experimental arrays are adequately described by a 
frustrated XY model, as assumed here. Of course, in the calculations, we are unable to vary 
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/ continuously as is possible in the experimental geometry. Thus, we are able to identify 
"peaks" only by considering nearby values of / not satisfying the relation / = |(1 — 1/N), 
and showing that these have smaller values of 7 than the / values in the chosen series. We 
have shown a few such values in Fig. 4. As one illustration, the value of 7 at / = 3/8 are 
larger than those at / = 17/48 and / = 19/48. 

Theron et al (JT2] account for the peaks at / = |(1 — 1/N) on the basis of a superlattice 
of vacancies in a background of the well-known checkerboard vortex ground state which 
forms at / = |. In their model, a change in the nature of the ground state is predicted at a 
critical frustration f c rs 0.468. When f c < f < |, a 2D vacancy superlattice is found to be 
energetically favored, while for 1/3 < / < f c , a vortex superlattice consisting of ID parallel 
stripes is preferred. 

In order to check this conjecture in our simulation, we have calculated the low- 
temperature vortex configuration of the vortex lattice at various values of /. In a Josephson 
junction array with frustration /, the vortex number n of each plaquette is defined as: 

£ (9i - 9j - An) + 2nf = 2un 

plaquette 

where the gauge invariant phase difference for each junction (#j — 9j — Ay) is in the range 
(— 7r, 7r] and the summation is taken in the counterclockwise direction. 

Some of our calculated vortex configurations are shown in Fig. 4. Fig. 4 (a) shows 
the configuration at / = 3/8 = 0.375 at temperature T = 0.01Ej/k B . It is easy to see 
there is indeed a striped domain structure in the vortex lattice, although the domain wall 



configuration is somewhat different from that predicted by ][L2|| . For lattices such that / > f c , 
it is necessary to carry out the calculations for relatively large lattices in order to distinguish 
the results from the / = 1/2 case. It proves very difficult to observe the predicted vacancy 
lattice structure. Fig. 4 (b) shows a representative low-temperature vortex configuration 
for an / value in this range, namely / = 17/36 = 0.472 at a temperature T = O.OlEj/ks- 
It is difficult to categorize the vortex lattice into any of the domain patterns suggested by 



12]. Nevertheless, we do obtain nearly all the peaks found experimentally, at most of the 



fields studied. We conclude that, while the experimental peaks do correspond to particularly 
stable vortex structures, there are many structures at these frustrations which are nearly 
equally stable, and which could all give rise to the observed peaks. 

Finally, we have checked for hysteresis at / = 1/6 in order to test whether the melting 
transition at that field is first-order. In 18 x 18 triangular arrays, we have used a Monte Carlo 
annealing method first to cool the array from a temperature l.bEj/ks to O.lEj/kB gradually, 
and then to reheat it. We find a clear hysteresis in the energy in the temperature range 
between 0.25Ej/kB and 0.32Ej/kB [cf. Fig. 5], suggestive of a first-order transition. This 
result is consistent with expectations based on three-dimensional calculations on stacked 
triangular lattices by Hetzel et al |fT7| , at / = 1/6. 
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V. CONCLUSION. 



In this paper, we have developed a new method to calculate the kinetic inductance 
of a Josephson junction array at zero temperature. This method maps a junction array 
onto a linear complex impedance network. The mapping permits us to calculate the zero- 
temperature kinetic inductance using much larger arrays than previously possible, provided 
that the ground state can be obtained by other means. The kinetic inductance calculated 
by this method agrees well with Monte Carlo calculations of the helicity modulus, in the 
low temperature limit. In both square and triangular arrays, the helicity modulus shows a 
series of peaks as a function of the frustration /. In triangular arrays, these peaks occur at 
/ = |(1 — 1/N), where N is an integer > 2, in agreement with experiment. 

The dynamical method discussed here has other potential applications. For example, 
it can be extended to treat arrays at finite frequencies, i. e., to study the variation of the 
effective admittance with frequency, simply by retaining the shunt resistances in Eq. (7). It 
can also be used with little change to treat underdamped arrays. Finally, it may be possible 
to extend this approach to finite temperatures, by including a Langevin noise term in the 
equations of motion for the phases. 
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Figure Captions 



1. Schematic diagrams of (a) an 8 x 8 square and (b) an 8 x 8 triangular Josephson junc- 
tion array. Each intersection represents a superconducting grain, which is connected 
to its nearest neighbors by Josephson coupling. In the dynamical inductance calcula- 
tion, the external current is injected into each grain at one boundary and extracted 
from the other side boundary; free boundary conditions are used in the direction of 
current injection, while periodic boundary conditions are used in the transverse di- 
rection. In the Monte Carlo calculations, periodic boundary conditions are applied in 
both directions. 

2. Calculated helicity modulus versus frustration / at four different temperatures in 
square arrays. Lines merely serve to guide the eye. For clarity, the curves for 
T = 0, 0.01Ej/kB and 0.05Ej/kB are shifted vertically by 0.3, 0.2, and 0.1 units 
respectively. The zero temperature result is obtained from a dynamical calculation of 
the inverse array inductance, using array sizes ranging from 120 x 120 to 180 x 180. 
The finite temperature results are obtained by Monte Carlo methods in 50 x 50 arrays 
for / = 1/5, 2/5, and in 48 x 48 arrays at all other values of / shown. 

3. Calculated helicity modulus versus frustration / at four different temperatures in tri- 
angular arrays, as obtained by Monte Carlo simulations. Lines merely serve to guide 
the eye. For clarity, we have vertically shifted the plots for T = O.OlEj/ks, O-ObEj/ks 
and 0.1Ej/k B by by 0.3, 0.2, and 0.1 units respectively. The array sizes are 56 x 56 
for / = 3/7; 50 x 50 for / = 1/5 and 2/5; and 48 x 48 for all other / values shown. 

4. Vortex configurations for triangular arrays as obtained by Monte Carlo annealing at 
temperature T = O.OlEj/ks at two different values of /: (a) / = 3/8, 32 x 32 array 
(the striped pattern can be seen by viewing the array from lower right to upper left); 
(b) / = 17/36, 36 x 36 array. Dark and light plaquettes denote the presence and 
absence of vortices. 

5. Energy per grain in 18 x 18 triangular arrays at / = 1/6, as calculated on cooling 
(solid disks) and on heating (open disks). The displayed results represent averages 
over five independent runs; lines serve to guide the eye. 
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Fig. 1 
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